///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file pfaffian.h
///
/// \brief   header file for functions related to pfaffian calculations for the generalized
///          pairing wavefunction
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_PFAFFIAN_HEADER
#define PCPS_PFAFFIAN_HEADER

#include <src/pcps.h>
#include <src/mantissa_exponent.h>

namespace pcps {

  /// \brief   Namespace for functions dealing with pfaffians and the pfaffian reference.
  namespace pfn {

    void init_occ_sites(const int np,
                        const int ns,
                        const int * const config,
                        int * const occ_sites,
                        int * const row_indices);

    template <class SCALAR> void build_occ_par_mat(const int np,
                                                   const int ns,
                                                   const SCALAR * const pairing_mat,
                                                   const int * const occ_sites,
                                                   SCALAR * const mat);

    pcps::MantExp<double> compute_pfaffian(const int n, double * const A, int * const iwork);

    pcps::MantExp<std::complex<double> > compute_pfaffian(const int n,
                                                          std::complex<double> * const A,
                                                          int * const iwork);

    template<class SCALAR> pcps::MantExp<SCALAR> get_pfaffian_value(const int np,
                                                                    const int ns,
                                                                    const int * const config,
                                                                    const SCALAR * const pairing_mat,
                                                                    int * const occ_sites,
                                                                    int * const iwork,
                                                                    SCALAR * const occ_par_mat);

    template <class SCALAR> void init_matrices(const int np,
                                               const int ns,
                                               const int * const config,
                                               const SCALAR * const pairing_mat,
                                               int * const occ_sites,
                                               int * const row_indices,
                                               SCALAR * const inv_mat,
                                               SCALAR * const work,
                                               int * const ipiv,
                                               int * const inv_perm);

    template <class SCALAR> SCALAR pfn_ratio(const int np,
                                             const int ns,
                                             const int nd,
                                             const int * const diff_sites,
                                             const int * const new_config,
                                             const int * const old_config,
                                             const SCALAR * const pairing_mat,
                                             SCALAR * const inv_mat,
                                             int * const occ_sites,
                                             int * const row_indices,
                                             SCALAR * const wv1);

    template <class SCALAR> void update(const int np,
                                        const int ns,
                                        const int nd,
                                        const int * const diff_sites,
                                        const int * const new_config,
                                        const int * const old_config,
                                        const SCALAR * const pairing_mat,
                                        SCALAR * const inv_mat,
                                        int * const occ_sites,
                                        int * const row_indices,
                                        SCALAR * const wv1,
                                        SCALAR * const wv2,
                                        SCALAR * const wv3);

    template <class SCALAR> void der_ratios(const int np,
                                            const int ns,
                                            const SCALAR * const inv_mat,
                                            const int * const occ_sites,
                                            const int * const row_indices,
                                            SCALAR * const ratios);

    template <class SCALAR> void make_first_row_col(const int np,
                                                    const int ns,
                                                    const int r,
                                                    SCALAR * const inv_mat,
                                                    int * const occ_sites,
                                                    int * const row_indices,
                                                    SCALAR * const wv1);

    void get_diff_occs(const int np,
                       const int ns,
                       const int nd,
                       const int * const diff_sites,
                       const int * const new_config,
                       const int * const old_config,
                       int * const new_occ_sites,
                       int * const old_occ_sites);

  } // end namespace pfn

} // end namespace pcps

#endif
